#Replication of the analysis in Hrbková, Kudrnáč (2024) "Fear, trust, 
# and compliance with covid-19 measures: A study of the mediating effect of 
#trust in government on the relationship between fear and compliance"

library(haven)
library(mediation)
library(dplyr)
w1 <- read_sav("/Users/lenkahrbkova/Documents/Papery 2023/JPP/wave1.sav")
w2 <- read_sav("/Users/lenkahrbkova/Documents/Papery 2023/JPP/wave2.sav")

######### WAVE 1 ###

# Create a new dataset with complete cases for all variables used in both models
w1_complete <- na.omit(w1[, c("comply_allNEW", "fear1", "fear3", "fear6", "trust_vlada", "AgeRec", "Gender", "Edu", "sympANO")])
w1_complete <- w1_complete %>%
  mutate(across(c(Gender, Edu), ~haven::as_factor(.))) %>%
  mutate(across(where(is.numeric), ~as.numeric(as.character(.))))  # Convert all other numeric columns to numeric
levels(w1_complete$Edu)

# Recoding the 'Edu' variable to binary
w1_complete$Edu <- ifelse(w1_complete$Edu == "Vysokoškolské", 1, 0)
table(w1_complete$Edu)

# Model 1 - Linear regression: fears - compliance
model1 <- lm(comply_allNEW ~ fear1 + fear3 + fear6 + trust_vlada + AgeRec + Gender + Edu + sympANO, data = w1_complete)
summary(model1)

# mediation models with the fear1
med_model1 <- lm(trust_vlada ~ fear1 + AgeRec + Gender + Edu + sympANO, data = w1_complete)
out_model1 <- lm(comply_allNEW ~ fear1 + trust_vlada + AgeRec + Gender + Edu + sympANO, data = w1_complete)
med_out1 <- mediate(med_model1, out_model1, treat = "fear1", mediator = "trust_vlada", sims = 1000)
summary(med_model1)
summary(out_model1)
summary(med_out1)

# models with the fear3
med_model2 <- lm(trust_vlada ~ fear3 + AgeRec + Gender + Edu + sympANO, data = w1_complete)
out_model2 <- lm(comply_allNEW ~ fear3 + trust_vlada + AgeRec + Gender + Edu + sympANO, data = w1_complete)
med_out2 <- mediate(med_model2, out_model2, treat = "fear3", mediator = "trust_vlada", sims = 1000)
summary(med_model2)
summary(out_model2)
summary(med_out2)

#models with fear6
med_model3 <- lm(trust_vlada ~ fear6 + AgeRec + Gender + Edu + sympANO, data = w1_complete)
out_model3 <- lm(comply_allNEW ~ fear6 + trust_vlada + AgeRec + Gender + Edu + sympANO, data = w1_complete)
med_out3 <- mediate(med_model3, out_model3, treat = "fear6", mediator = "trust_vlada", sims = 1000)
summary(med_model3)
summary(out_model3)
summary(med_out3)


## test of normality WAVE 1
# Extract residuals from the mediator model
med_model1_resid <- resid(med_model1)
med_model2_resid <- resid(med_model2)
med_model3_resid <- resid(med_model3)

# Perform Shapiro-Wilk test for normality on these residuals
shapiro_test_med_model1 <- shapiro.test(med_model1_resid)
shapiro_test_med_model2 <- shapiro.test(med_model2_resid)
shapiro_test_med_model3 <- shapiro.test(med_model3_resid)
# Print the result
print(shapiro_test_med_model1)
print(shapiro_test_med_model2)
print(shapiro_test_med_model3)

# Extract residuals from the outcome model
out_model1_resid <- resid(out_model1)
out_model2_resid <- resid(out_model2)
out_model3_resid <- resid(out_model3)

# Perform Shapiro-Wilk test for normality on these residuals
shapiro_test_out_model1 <- shapiro.test(out_model1_resid)
shapiro_test_out_model2 <- shapiro.test(out_model2_resid)
shapiro_test_out_model3 <- shapiro.test(out_model3_resid)

# Print the result
print(shapiro_test_out_model1)
print(shapiro_test_out_model2)
print(shapiro_test_out_model3)

###bootstrapping
# Run mediation with bootstrapping
med_out1boot <- mediate(med_model1, out_model1, treat = "fear1", mediator = "trust_vlada", sims = 1000, boot = TRUE)
med_out2boot <- mediate(med_model2, out_model2, treat = "fear3", mediator = "trust_vlada", sims = 1000, boot = TRUE)
med_out3boot <- mediate(med_model3, out_model3, treat = "fear6", mediator = "trust_vlada", sims = 1000, boot = TRUE)

# Summary of the mediation analysis with bootstrapped confidence intervals
summary(med_out1boot)
summary(med_out2boot)
summary(med_out3boot)


############ WAVE 2 #######################
w2_complete <- na.omit(w2[, c("comply_allNEW", "fear1", "fear3", "fear6", "trust_vlada", "AgeRec", "Gender", "Edu", "sympANO")])
w2_complete <- w2_complete %>%
  mutate(across(c(Gender, Edu), ~haven::as_factor(.))) %>%
  mutate(across(where(is.numeric), ~as.numeric(as.character(.))))  # Convert all other numeric columns to numeric

levels(w2_complete$Edu)
# Inspect the unique values and types in the original Edu variable
unique(w2$Edu)
# Convert Edu to a factor using the appropriate method
w2_complete <- w2 %>%
  mutate(Gender = as.factor(Gender),
         Edu = factor(Edu))  

w2_complete$Edu <- ifelse(w2_complete$Edu == "4", 1, 0)
table(w2_complete$Edu)
levels(w2_complete$Edu)

# Model 2 - Linear regression: fears - compliance
model2<- lm(comply_allNEW ~ fear1 + fear3 + fear6 + trust_vlada + AgeRec + Gender + Edu + sympANO, data = w2_complete)
summary(model2)

##### models with fear1
med_model1_w2 <- lm(trust_vlada ~ fear1 + AgeRec + Gender + Edu + sympANO, data = w2_complete)
out_model1_w2 <- lm(comply_allNEW ~ fear1 + trust_vlada + AgeRec + Gender + Edu + sympANO, data = w2_complete)
med_out1_w2 <- mediate(med_model1_w2, out_model1_w2, treat = "fear1", mediator = "trust_vlada", sims = 1000)
summary(med_model1_w2)
summary(out_model1_w2)
summary(med_out1_w2)

# models with the fear3
med_model2_w2 <- lm(trust_vlada ~ fear3 + AgeRec + Gender + Edu + sympANO, data = w2_complete)
out_model2_w2 <- lm(comply_allNEW ~ fear3 + trust_vlada + AgeRec + Gender + Edu + sympANO, data = w2_complete)
med_out2_w2 <- mediate(med_model2_w2, out_model2_w2, treat = "fear3", mediator = "trust_vlada", sims = 1000)
summary(med_model2_w2)
summary(out_model2_w2)
summary(med_out2_w2)

#models with fear6
med_model3_w2 <- lm(trust_vlada ~ fear6 + AgeRec + Gender + Edu + sympANO, data = w2_complete)
out_model3_w2 <- lm(comply_allNEW ~ fear6 + trust_vlada + AgeRec + Gender + Edu + sympANO, data = w2_complete)
med_out3_w2 <- mediate(med_model3_w2, out_model3_w2, treat = "fear6", mediator = "trust_vlada", sims = 1000)
summary(med_model3_w2)
summary(out_model3_w2)
summary(med_out3_w2)

## test of normality in WAVE 2
# Extract residuals from the mediator model
med_model1_w2_resid <- resid(med_model1_w2)
med_model2_w2_resid <- resid(med_model2_w2)
med_model3_w2_resid <- resid(med_model3_w2)

# Perform Shapiro-Wilk test for normality on these residuals
shapiro_test_med_model1_w2 <- shapiro.test(med_model1_w2_resid)
shapiro_test_med_model2_w2 <- shapiro.test(med_model2_w2_resid)
shapiro_test_med_model3_w2 <- shapiro.test(med_model3_w2_resid)
# Print the result
print(shapiro_test_med_model1_w2)
print(shapiro_test_med_model2_w2)
print(shapiro_test_med_model3_w2)

# Extract residuals from the outcome model
out_model1_w2_resid <- resid(out_model1_w2)
out_model2_w2_resid <- resid(out_model2_w2)
out_model3_w2_resid <- resid(out_model3_w2)

# Perform Shapiro-Wilk test for normality on these residuals
shapiro_test_out_model1_w2 <- shapiro.test(out_model1_w2_resid)
shapiro_test_out_model2_w2 <- shapiro.test(out_model2_w2_resid)
shapiro_test_out_model3_w2 <- shapiro.test(out_model3_w2_resid)

# Print the result
print(shapiro_test_out_model1_w2)
print(shapiro_test_out_model2_w2)
print(shapiro_test_out_model3_w2)


## Bootstrap mediation models
med_out1boot_w2 <- mediate(med_model1_w2, out_model1_w2, treat = "fear1", mediator = "trust_vlada", sims = 1000, boot = TRUE)
med_out2boot_w2 <- mediate(med_model2_w2, out_model2_w2, treat = "fear3", mediator = "trust_vlada", sims = 1000, boot = TRUE)
med_out3boot_w2 <- mediate(med_model3_w2, out_model3_w2, treat = "fear6", mediator = "trust_vlada", sims = 1000, boot = TRUE)

# Summary of the mediation analysis with bootstrapped confidence intervals
summary(med_out1boot_w2)
summary(med_out2boot_w2)
summary(med_out3boot_w2)


## Correction for multiple testing
p_values1 <- c(2.09e-12, 0.605631, 0.017972, 1.29e-06, 6.96e-07, 9.75e-07, 0.000506, 0.393415,
               1.28e-07, 0.0896, 0.0594, 0.3242, 1e-16, 2.88e-16, 1.67e-07,
               2.19e-07, 1.82e-06, 0.00061, 0.64834, 1e-16, 1e-16,
               1e-16, 1e-16, 0.09416, 0.00538, 0.16921, 0.59228,
               1e-16, 7.15e-06, 2.33e-10, 9.60e-12, 1.45e-07, 0.00437, 0.86415,
               0.12, 1e-16, 1e-16, 0.12, 0.00302, 0.01242,
               0.36971, 0.78169, 1e-16, 0.8940, 7.94e-11, 1.98e-10, 1.32e-08,
               0.0189, 0.9944, 0.01, 0.90, 0.60, 0.59)

p_values2 <- c(2e-16, 0.0565, 1.68e-12, 2.69e-06, 2.31e-08, 3.09e-11, 
               0.5000, 0.011, 2.65e-06, 0.00874, 0.03982, 0.00488, 2e-16, 
               2e-16, 4.96e-12, 9.50e-08, 9.33e-10, 0.8677, 0.0166, 
               2e-16, 2e-16, 2e-16, 2e-16, 0.00297, 0.00455, 0.49804, 
               0.01764, 2e-16, 0.00146, 2e-16, 3.09e-12, 7.02e-13, 0.78863,
               0.12623, 0.002, 0.004, 0.038, 0.040, 2e-16, 0.00265, 0.92628,
               0.11771, 2e-16, 5.85e-10, 3.90e-09, 5.47e-11, 2e-16, 0.259, 
               0.164, 2e-16, 2e-16, 2e-16, 2e-16)

# Applying BY correction
p_adjusted_BY <- p.adjust(p_values1, method = "BY")
p_adjusted_BY2 <- p.adjust(p_values2, method = "BY")#wave2

# Rounding the adjusted p-values to three decimal places
p_adjusted_by_rounded <- round(p_adjusted_BY, 3)
p_adjusted_by2_rounded <- round(p_adjusted_BY2, 3)#wave2

# Replacing very small p-values with "<2e-16" in the final output
p_adjusted_by_final <- ifelse(p_adjusted_by_rounded < 2e-16, "<2e-16", p_adjusted_by_rounded)
p_adjusted_by2_final <- ifelse(p_adjusted_by2_rounded < 2e-16, "<2e-16", p_adjusted_by2_rounded)

# Print the adjusted and rounded p-values
print(p_adjusted_by_final)
print(p_adjusted_by2_final)

